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1 ABSTRACT 

EUPDF-II provides the solution for the species 
and temperature fields based on an evolution equa- 
tion for PDF (Probability Density Function) and it 
is developed mainly for application with sprays, com- 
bustion, parallel computing and unstructured grids. 
It is designed to be massively parallel and could eas- 
ily be coupled with any existing gas-phase CFD and 
spray solvers. The solver accommodates the use of an 
unstructured mesh with mixed elements of either tri- 
angular, quadrilateral, and/or tetrahedral type. The 
manual provides the user with an understanding of 
the various models involved in the PDF formulation, 
its code structure and solution algorithm, and various 
other issues related to parallelization and its coupling 
with other solvers. The source code of EUPDF-II will 
be available with National Combustion Code (NCC) 
as a complete package. 

2 NOMENCLATURE 

A pre-exponent coefficient in an Arrhenius 
reaction rate term 

a non-unity exponent in an Arrh- 
enius reaction rate term 
a n outward area normal vector 
of the nth surface, m 2 
Bk Spalding transfer number 
b non-unity exponent in an 

Arrhenius reaction rate term 
C p specific heat, J / (kg K) 

C<$> the constant used in Eq. (7) 
c n convection /diffusion coefficient 
of the nth face, kg/s 

D turbulent diffusion coefficient, m 2 /s 
E a activation energy in an Arrhenius 
reaction rate term, J/kg-mole 
h specific enthalpy, J/kg 


Jf diffusive mass flux vector, kg/ms 

k turbulence kinetic energy, m 2 /s 2 

Ik mixture latent heat of evaporation, J/kg 
Ik.eff effective latent heat of evaporation, 

J/kg (defined in Eq. (11)) 
l kin heat of vaporization at normal 

boiling point, J/kg 
Mi molecular weight of ith 

species, kg/kg-mole 
rrik droplet vaporization rate, kg/s 

N av number of time steps employed in 

the PDF time-averaging scheme 
Nf number of surfaces contained in 

a given computational cell 
N m total number of Monte Carlo 

particles per grid cell 

N p total number of computational cells 

rik number of droplets in kth group 

P pressure, N/m 2 

P c critical pressure, N/m 2 

P r Prandtl number 

p joint scalar PDF 

R u universal gas constant, 8314.4 J/(kg-mole K) 

Re Reynolds number 

Smic liquid source contribution of the 

gas-phase continuity equation 
Smie liquid source contribution of the 

gas-phase energy equation 
s m i m liquid source contribution of the 

gas-phase momentum equations 
Sfnis liquid source contribution of the 

gas-phase species equations 

s a liquid source contribution for the a-th variable 

T temperature, K 

T nb normal boiling point, K 

T c critical temperature, K 
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, partial differentiation with respect 
to the variable followed by it 
a the total number of composition 
variables in Eq. (5) 

Superscripts 


t time, s 

m velocity, m/s 

V c volume of the computational cell, m 3 
w a represents the chemical reaction rate term 
in Eq. (5), 1/s 

Wj represents the chemical reaction rate term 
in Eq. (3), 1/s 

Cartesian coordinate in the ith direction, m 
yj mass fraction of jth species 

x spatial vector 

X mole fraction 

At local time step used in the 

PDF computations, s 

Atf local time step used in the flow solver, s 
AV computational cell volume, m 3 

5 Dirac-delta function 

e rate of turbulence dissipation, m 2 /s 3 

Cj fractional mass evaporating rate of species 

at the droplet surface 

e as species mass fraction at the droplet surface 
r<£ turbulent diffusion coefficient, kg/ms 

A thermal conductivity, J/(ms K) 

fx dynamic viscosity, kg/ms 

cj turbulence frequency, 1/s 

(j> represents the solution field of the joint 
PDF transport equation 
) independent composition space 

p density, kg/m 3 

pi n liquid density (at 1 bar, 273.15 K), kg/m 3 
a dimensionality of ^0-space 

r stress tensor term, kg/ms 2 

6 void fraction 

Subscripts 

/ represents conditions associated with fuel 
g global or gas-phase 

i index for the spatial coordinate 

or the species component 
j index for the species component 

k droplet group or liquid-phase 

l liquid-phase 

m conditions associated with N m 

n nth-face of the computational cell 

o either initial conditions or oxidizer 

p conditions associated with the 

properties of a grid cell 
s represents conditions at the droplet 
surface or adjacent computational cell 
t conditions associated with time 

a index for the scalar component of 

the joint PDF equation 


Favre averaging 

time averaging or average based on the Monte 
Carlo particles present in a given cell 

* the solution step associated with Eq. (14) 

* * the solution step associated with Eq. (15) 

* * ★ the solution step associated with Eq. (16) 

// fluctuations 

3 INTRODUCTION 

The gas-turbine combustor flows are often char- 
acterized by a complex interaction between various 
physical processes such as the interaction between the 
liquid and gas phases, heat release associated with 
chemical kinetics, and radiative heat transfer asso- 
ciated with highly absorbing and radiating species. 
The rate controlling processes often interact with 
each other at various disparate time and length scales. 

In particular, turbulence plays an important role in 
determining the rates of mass transfer, heat transfer, 
chemical reactions, and droplet evaporation. Most of 
the turbulence closure models for reactive flows have 
difficulty in treating nonlinear reaction rates [1-2]. 
The use of assumed shape PDF methods was found to 
provide reasonable predictions for pattern factors and 
NOx emissions at the combustor exit. However, their 
extension to multi-scalar chemistry becomes quite in- 
tractable. The solution procedure based on the mod- 
eled joint scalar PDF transport equation has an ad- 
vantage in that it treats the nonlinear reaction rates 
without any approximation. In this approach, the ve- 
locity and turbulence fields are solved using a conven- 
tional CFD solver, a modeled PDF transport equa- 
tion provides the solution for the species and temper- 
ature fields, and an appropriate spray formulation is 
used for the liquid-phase representation. The method 
based on the combined PDF/CFD /spray calculations 
holds the promise of providing an accurate represen- 
tation for several important combustion phenomena 
associated with the prediction of emissions (CO and 
NOx), UHC (Unburnt HydroCarbons), and flame ex- 
tinction and blow-off limits [1]. 

However, one of the major disadvantages of 
this method is that all of the composition variables 
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also appear as the independent variables of the PDF 
transport equation. Because of the large dimension- 
ality, it is almost impossible to obtain a solution for 
the joint scalar PDF transport equation based on con- 
ventional finite-difference techniques [1]. However, 
Monte Carlo methods are better suited for obtaining 
a solution to this PDF equation because it was shown 
that the required computational effort increases only 
linearly with the increase in the dimensionality of the 
PDF transport equation [1], Other than some appli- 
cation to simple flows, there is a very limited evidence 
of the application of this method to the calculation 
of practical combustion flows [3]. It is because the 
Monte Carlo methods tend to be both computation- 
ally very time consuming and require a large com- 
puter memory for application to 3D flows [3]. How- 
ever, the success of any numerical tools used in multi- 
dimensional combustor modeling depends not only on 
the modeling- and numerical-accuracy requirements 
but also on the computational-efficiency considera- 
tions as determined by the computer memory and 
turnaround times afforded by the present-day com- 
puters. With the aim of demonstrating its viability 
to flows representative of a gas-turbine combustor, we 
have undertaken the task of extending this technique 
in a number of significant ways: 

1. In order to demonstrate the importance of chem- 
istry/turbulence interactions in the modeling 
of reacting sprays, we have extended the joint 
scalar Monte Carlo PDF approach to the mod- 
eling of spray flames. 

2. It facilitates the use of both unstructured grids 
and parallel computing and, thereby, facilitat- 
ing large-scale combustor computations involv- 
ing complex geometrical configurations. 

3. We have developed and implemented several 
numerical convergence techniques such as local 
time-stepping and various other averaging meth- 
ods in order to accelerate convergence to a steady 
state. 

4. The PDF method provided favorable agreement 
when applied to several supersonic diffusion 
flames as well as several other subsonic spray 
flames [4-9]. 

Some of our work on the Monte Carlo PDF 
method could be found in Refs. [4-9]. Initially, the 
emphasis of our work was on extending this tech- 
nique to the modeling of compressible reacting flows 


[4]. The Monte Carlo solver was used in conjunction 
with several density-based CFD codes for the mean- 
velocity and turbulence fields. Several averaging pro- 
cedures introduced in Refs. [4-5] proved to be use- 
ful in providing smooth Monte Carlo solutions to the 
CFD solver. The PDF method provided favorable 
results when applied to several supersonic diffusion 
flames [4-5]. 

Later on this approach was further extended to 
the modeling of spray flames and parallel comput- 
ing and, thereby, combining the novelty of the PDF 
method with the ability to run on parallel architec- 
tures [6]. This algorithm was implemented on the 
Cray T3D at NASA Lewis Research Center, a mas- 
sively parallel computer, with an aggregate of 64 Pro- 
cessor Elements (PEs). The computer code was writ- 
ten in Cray MPP (Massively Parallel Processing) For- 
tran. The application of this method to both open as 
well as confined axisymmetric swirl-stabilized spray 
flames showed reasonable agreement with the mea- 
sured data [6]. 

With the aim of advancing the current multi- 
dimensional computational tools used in the design 
of advanced technology combustors, two new com- 
puter codes, LSPRAY - a Lagrangian spray solver 
[10] and EUPDF - an Eulerian Monte Carlo PDF 
solver [11], were developed, thereby extending our 
previous work [6] on the Monte Carlo PDF and sprays 
to unstructured grids as a part of NCC development. 
The combined unstructured 3D CFD / S pray / Mont e- 
Carlo-PDF solver is designed to be massively par- 
allel and accommodates the use of an unstructured 
mesh with mixed elements comprised of either trian- 
gular, quadrilateral, and/or tetrahedral type. This 
was done mainly to facilitate representation of com- 
plex geometries with relative ease [7-9]. Also, it is 
noteworthy that considerable effort usually goes into 
the grid generation of practical combustion devices 
which tend to be very complex in both shape and 
configuration. The grid generation time could be re- 
duced considerably by making use of existing unstruc- 
tured grid generators such as GRIDGEN. A current 
status of the use of the parallel computing in tur- 
bulent reacting flows involving sprays, scalar Monte 
Carlo PDF and unstructured grids was described in 
Ref. [7]. It outlined several numerical techniques de- 
veloped for overcoming some of the high computer 
time-and-storage requirements associated with the 
use of Monte Carlo solution methods. The parallel 
performance of both the PDF and CFD computa- 
tions was found to be excellent but the results on the 
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spray computations showed reasonable performance. 
Our preliminary estimates indicated that it was well 
within the reach of modern parallel computer’s ca- 
pacity to perform a realistic gas-turbine combustor 
simulation [6-9]. 

Some salient features of the Monte Carlo PDF 
module are summarized below: 

1. The scalar Monte Carlo PDF solver has been 
extended to mixed unstructured grid elements 
of triangular, quadrilateral, tetrahedron, wedge, 
and hexahedron elements. 

2. EUPDF-II is currently coupled with an unstruc- 
tured flow (CFD) solver of NCC [12-13], and a 
Lagrangian spray solver - LSPRAY [10], which 
were selected to be as the integral components 
of the NCC cluster of modules. 

3. The PDF code receives the mean velocity and 
turbulence fields from the flow solver. And, 
also, it receives the source terms arising from the 
liquid-phase contribution from the spray solver if 
needed. 

4. The PDF solver provides the species and tem- 
perature solution to the CFD solver and, also, 
to the spray solver if needed. 

5. The PDF code is written in Fortran 77/90 with 
PVM or MPI calls for parallel computing. It 
enables the computations to be performed with 
equal ease on both massively parallel computers 
as well as workstation clusters. 

The furnished code demonstrates the success- 
ful methods used for parallelization and coupling of 
the PDF to the flow code. First, complete details 
of the Monte Carlo PDF solution procedure is pre- 
sented along with several other numerical issues re- 
lated to the coupling between the CFD, EUPDF-II, 
and LSPRAY modules. It is followed by a brief de- 
scription of the combined parallel performance of the 
three solvers (CFD, EUPDF-II, and LSPRAY) along 
with a brief summary of the validation cases. 

4 GOVERNING EQUATIONS FOR THE 
GAS PHASE 

Here, we summarize the conservation equations 
for the gas-phase in Eulerian coordinates [14]. These 
equations are valid for a dilute spray with a void frac- 
tion of the gas, 0, close to unity. The void fraction is 


defined as the ratio of the equivalent volume of gas 
to a given volume of a gas and liquid mixture. This 
is done for the purpose of identifying the interphase 
source terms arising from the exchanges of mass, mo- 
mentum, and energy with the liquid-phase. 

The conservation of the mass leads to: 

[pVc],t T [pLcU-j]^ — Smlc “ ^ ^ ^ k (I) 

k 

For mass conservation, the source term is given as 
a summation over different classes of droplets. Each 
class represents the average properties of a different 
polydisperse group of droplets. Here, n* represents 
the number of droplets in a given class and ra^ rep- 
resents the corresponding mass vaporization rate. 

For the conservation of the jth species, we have: 

[p^c2/j],t 4“ [pVc'U‘iyj] i xi [pUc^2/jf,«t] 

—pV c ibj = s mls = ^ ej n k m k (2) 

k 

where 

Wj — 0 and ^ tj — 1 

3 3 

For the species conservation, the source term con- 
tains an additional variable, ej, which is defined as 
the fractional vaporization rate for species j. 

For the momentum conservation, we have: 

[pV c Uj], t + [pV c UiUj} tXj + \pV c } )Xi ~. 

— [(1 — @)Vc.1~Uj] ,xj = Smlrrt = 

n k m k u ki -Y, y Pk r l n k umj (3) 

k k 

where the shear stress in Eq. (3) is given by: 

Tij — + Uj iXi ] — 

For the momentum conservation, the first source 
term represents the momentum associated with liquid 
fuel vapor and the second represents the momentum 
change associated with droplet drag. 

For the energy conservation, we have: 
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[^ e A],*+[^ e u i A],« 4 -[tfv e Ar, ie< ] 1 « < -[(i-tf)v c A/r 1 , < ] lie4 

-[QV e p\,t = Sm< e = 'Y^n k m k ( h s - h,e}}) (4) 

k 

Similarly, the energy conservation has a source first 
term associated with liquid fuel vapor and the second 
represents the heat loss associated with the latent 
heat of vaporization and it also contains an additional 
heat flux (loss or gain) to the droplet interior from the 
ambient. 

The main purpose of the spray solver is to cal- 
culate the source terms arising from the exchanges 
of the mass, momentum, and energy and, then, feed 
that appropriate information to the gas-phase solver. 
In the case of NCC, it supplies the source terms to 
the CFD solver and, also, to the Monte Carlo PDF 
solver if needed. 

5 GAS-PHASE SCALAR JOINT PDF 
EQUATION 

The transport equation for the density- weighted 
joint PDF of the compositions, p, is: 

p\p],t + pUi[P}, Xi 
{Transient} {Mean convection} 

+ [pWai^p],^ = -\p < u" I rp_ > p\, Xi 
{Chemical reactions} {Turbulent convection} 

-[? < -Jtxi 1 1 > pUc ~\P< \ s <* 1 1 > ( 5 ) 

r r 


{Molecular mixing} 

{Liquid— phase contribution} 

where 

w a 

= chemical source term 
for the a-th 
composition variable, 

< u” I Ip_ > 

= conditional average of 
Favre velocity 
fluctuations, 

< 7 J txi\±> 

= conditional average 
of scalar dissipation, & 

< j*a \±> 

— conditional average of 
spray source terms. 


The terms on the left-hand side of the above 
equation could be evaluated without any approxima- 
tion, but the terms on the right-hand side of the equa- 
tion require modeling. The first term on the right rep- 
resents transport in physical space due to turbulent 


convection [1]. Since the joint PDF, p, contains no in- 
formation on velocity, the conditional expectation of 
< u” | ip > needs to be modeled. It is modeled based 
on a gradient-diffusion model with information sup- 
plied on the turbulent flow field from the flow solver 

[i]- 

- < u" I tp>p = Y^p tXi (6) 

The fact that the turbulent convection is modeled 
as a gradient-diffusion makes the turbulent model no 
better than the k—e model. However, it is noteworthy 
that the option of using a non-linear k — e turbulence 
model of Shih et al. [15], seems to alleviate some of 
the modeling uncertainties associated with the use of 
the standard k — e model. 

The second term on the right-hand side repre- 
sents transport in the scalar space due to molecu- 
lar mixing. A mathematical description of the mix- 
ing process is rather complicated, and the interested 
reader is referred to Ref. 2. Molecular mixing is ac- 
counted for by making use of the relaxation to the 
ensemble mean submodel [2]. 

< -• J t Xi 1 1 >= - <A«) (7) 

where u — e/&, and C<t> is a constant. For a 

conserved scalar in a homogeneous turbulence, this 
model preserves the PDF shape during its decay, but 
there is no relaxation to a Gaussian distribution [1]. 
However, the results of Ref. 16 indicate that the 
choice between different widely-used mixing models 
is not critical in the distributed reaction regime of 
premixed combustion as long as the turbulent mix- 
ing frequencies are above 1000 Hz. Most practical 
combustors seem to operate at in-flame mixing fre- 
quencies of 1000 Hz and above. The application of 
this mixing model seemed to provide some satisfac- 
tory results when applied to flows representative of 
those encountered in the gas-turbine combustion [16]. 

The third term on the right-hand side represents 
the contribution from the spray source terms: 

< L a | V>> = -~y'T, n kTn k {e as - <f> a ) (8) 

where <j> a — y a when a = 1, 2, ..., a — 1 

< | ip >= + h ks - (pa) 

P P (9) 
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where <t> a =o — h and is defined by: 


where 


O — 1 

h = ^2 yihi 

%= 1 


(10) 


hi — 



C P i = Au + A 2i T + A 3i T 2 + A 4i T 3 + A 5i T% 
vv i 

h°^ i is the heat of formation of ith species, R u is 
the universal gas constant, e as is the fractional mass 
evaporating rate of species at the droplet surface, and 
l keff is the effective latent heat of vaporization as 
modified by the heat loss to the droplet interior: 


where the subscript n refers to the nth-face of the 
computational cell. The coefficient c n represents the 
combined flux in physical space due to convection and 
diffusion through the nth-face of the computational 
cell, p. The convection /diffusion coefficients in the 
above equation are determined by one of the following 
two expressions: 


. = r*( 


-a r 


A14> + A V s 


) + max[ 0, —pan-lLr, 


c n = maa;[jG.5po„.w n |, ) 1 ~ ®Apa n .u n 


and 


£' 


= h + 4 %jf(fr) < ( u > 

where the mixture latent heat of vaporization, /&, is 
given by 


In both the above expressions for c n , a cell-centered 
finite-volume derivative is used to describe the vis- 
cous fluxes; but an upwind differencing scheme is used 
for the convective fluxes in the first expression and a 
hybrid differencing scheme in the second. 


and 


Ik — ^ ^ 


^ki — ^kin ^ 


T ci — T \°- 38 
T ci -T bi ) 5 


6.1 Numerical Method Based on 
Approximate Factorization 

The transport equation is solved by making use of an 
approximate factorization scheme [1]. Eq. (12) can 
be recast as: 


rp IkinMj/ R u 

IfrinMi/ {R u tbni) — 

More details on Eq. (11) can be found in a 
companion report on LSPRAY-II: a Lagrangian spray 
module. 


Pp(±,t + At) = 

(I+AtR){I+At.S){I+AtM){I+AtT)p p (±,t)+0{At 2 ) 

(13) 


6 PDF SOLUTION ALGORITHM 

Partial finite-differencing and rearrangement of 
the terms associated with the convection and diffu- 
sion of Eq. (5) yields: 


where / represents the unity operator and T, M, S, 
and R denote the operators associated with spatial 
transport, molecular mixing, spray, and chemical re- 
actions, respectively. The operator is further split 
into a sequence of intermediate steps: 


c /\ f 

Pp(±, t + At) = (1 _ P AV )p P {±, t ) 

= (I + A tT)p p {rp,t) 

(14) 

+ ] -^yPni^A) ~ A t[w a (^)p] t ^ a 

n ^ 

p;*(±,t) = (I + AtM)p* p (±,t) 

(15) 

At[< -j& t \±>PU«~ A *[< y a \±> pu» 

P ’ (12) 

p;**(p L ,t) = {i + Ats)p;*(±,t) 

(16) 


NASA/CR— 2004-2 1 3073 


6 



P p (& t + At) = (I + A tR)p™{±, t) (17) 

The operator-splitting method is used to provide the 
solution for p by making use of a Monte Carlo tech- 
nique. In the Monte Carlo simulation, the continuos 
PDF is replaced by a a discrete PDF which is defined 
in terms of a large ensemble of stochastic particles. 
The ensemble-averaged PDF over N m delta functions 
replaces the average based on a continuous PDF [1]. 

i 

Ppm(4 0 =< P P {4 0 >= ttY ~ 0 ") ( 18 ) 

1 m n= 1 

The discrete PDF p pm (V0 is defined in terms of an 
ensemble of N m Dirac- delta functions each associated 
with a corresponding set of (f) n : n — 1, 2, 3 .. JV m . The 
statistical error in this approximation is proportional 
to iV m 1/2 . 

Using the operator-splitting method, the steps 
associated with the spatial transport due to convec- 
tion and turbulent diffusion as well as the steps asso- 
ciated with the transport in the scalar space due to 
molecular mixing, spray contribution, and chemical 
kinetics are advanced in a series of sequential steps 
as given by Eqs. (14)-(17). 

6.2 Convection/Diffusion Step 

The first step associated with convection /diffusion is 
given by: 


P* p {±, t) = {I + A tT)p p {ip, t ) = 

(! - V+Y *) ( 19 ) 

This step is simulated by replacing a randomly se- 
lected number of particles as given by <£*(£), ..,<^ (t) 
(Nt = the nearest integer of Cn ^y~ ) with by a set of 
the same number of randomly selected particles cho- 
sen from the neighboring cell located on the opposite 
side of the corresponding n-th face. 

6.3 Numerical Issues Associated With Fixed 
Versus Variable Time Step 

It is obvious from the above equation that a nec- 
essary criterion for stability requires satisfaction of 
< 1. When the computations are performed 
with a fixed time step, this criterion tends to be too 
restrictive for most applications. Depending on the 


flow configuration, the allowable maximum time in- 
crement At is likely to be limited by a region of the 
flow field where convective fluxes dominate (such as 
close to injection holes). But in the main stream, the 
flow is usually characterized by much lower velocities. 
Resolution considerations require a higher concentra- 
tion of the grid in certain regions of the flow r -field than 
in others. For example, more grid cells are clustered 
in regions where boundary layers are formed. In such 
regions the allowable maximum time increment might 
be limited in a direction dominated by the largest of 
the diffusive fluxes as determined by T^/Ax. This 
problem gets magnified if the cells also happen to be 
highly skewed. 

Such restrictions on the allowable maximum 
time step could lead to a frozen condition when the 
Monte Carlo simulation is performed with a limited 
number of stochastic particles per cell. For clarity, 
let us consider the following criterion: 


N m > 


pAV 

c n At 


( 20 ) 


which has to be satisfied at all grid nodes. It is es- 
timated that about 10 3 stochastic particles per cell 
are needed in order to avoid the so-called frozen con- 
dition for performing a typical 3-D gas-turbine com- 
bustor calculation. The frozen condition is referred to 
as a state in which no transfer of stochastic particles 
takes place between the neighboring cells when N m 
falls below a minimum required. Scheurlen et al. [3] 
were the first ones to recognize the limitations asso- 
ciated with the use of a fixed time step in the Monte 
Carlo PDF computations. 

However, our experience has shown that this 
problem can be overcome by introducing the con- 
cept of local time-stepping which is a convergence 
improvement technique widely used in many of the 
steady-state CFD computations. In this approach, 
the solution is advanced by making use of a local 
time step which could be different at each one of its 
computational nodes. In our present computations, 
it is determined based on 


At = min(CifAtf, 7^-) (21) 

where C t f and Ct are calibrated constants, Atf is the 
local time step obtained from the flow module. The 
local time step is chosen such that it permits trans- 
fer of enough particles across the boundaries of the 
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neighboring cells while ensuring that the time step 
used in the PDF computations does not deviate con- 
siderably from the time step used in the CFD solver. 

6.4 Molecular Mixing Step 

The second step associated with molecular mixing is 
given by: 

^ ( 22 ) 

where C# is an empirical constant with a value of 
about unity. 

The solution for this equation is updated by: 

<pz = <Pi + (r a -r a )e~ c *“ At ( 23 ) 

6.5 Spray Step 

The third step associated with the spray contribution 
is given by: 


Finally, the fourth step associated with chemical re- 
actions is given by: 

^ (28) 

at p W f \\ o 

where <j> a = 2//- 

^ = - Uo ^AA a (^) b e-^ (29) 

at p Wf Wo 

where <j> a = y 0 - 

= 0 (30) 

dt 

where 4> a = h. 

The numerical solution for Eqs. (28)-(30) is ob- 
tained by an implicit Euler scheme [17]. The result- 
ing non-linear algebraic equations are solved by the 
method of quasi-linearization [18]. 

6.7 Details of Combustion Chemistry 


d<t>a 

dt 


1 

pAV 


+ " $<*) 


where 4> a = y a when a = 1,2, ..., a — 1 


(24) 


= TXT7 y ~] n k‘>nk(—h,eff + hks ~ <t>a) (25) 
dt pAV ' 

where 4> a — h when a = a The solution for the above 
equations is upgraded by a simple explicit scheme: 


K 


At 52 n k m k 


pAV 


+ - 


A t52nkmk 


pAV 


) 


where a < a — 1 


(26) 


,*** AtJ2n k m k | u \ 

<P a = + h k3 ) 

. A t 52 n k m k , 

+ ^ (1 “ pAV 


In this section, we present an example of how com- 
bustion chemistry is handled for the case of n-heptane 
when it is modeled by a single-step global mecha- 
nism of Westbrook and Dryer [19]. The correspond- 
ing rate constants in Eqs. (28)-(29) are given by 
A = 0.286E + 10 , a = 0.25, b - 1.25, and 
E a — 0. 151^ + 05. This global combustion mecha- 
nism was reported to provide adequate representation 
for temperature histories in flows not dominated by 
long ignition delay times. For example, the overall re- 
action representing the oxidation involving n-heptane 
is given by: 


C 7 l?i 6 + 11(02 + 3.767V 2 )-+ 

7002 + SH 2 0 + 41.36 N 2 (31) 

Because of the constant-Schmidt-number as- 
sumption made in the PDF formulation and based on 
atomic balance of the constituent species, the mass 
fractions of JV 2 , 00 2 , and H 2 0 can be shown to be 
related to the mass fractions of 0 2 and 0 7 i7i6 by the 
following expressions: 


where a — a. After a new value for enthalpy is up- 
dated, the temperature is determined iteratively from 
the solution of Eq. (10). 

6.6 Reaction Step 


y H2 o -K 2 - K\ K 2 yo 2 ~ K 2 yc 7 H 16 

yco 2 — K3IJH2O (32) 

yjv 2 = 1 — K 2 - K 2 K z - yo 2 (l — K\K 2 - KiK 2 Kz)- 
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yC 7 H r6 {l — K 2 — K 2 K 3 ) 

where K x = 4.29, K 2 = 0.08943, and K 3 = 2.138. 

Making use of Eq. (32) leads to a considerable 
savings in computational time as it reduces the num- 
ber of variables in the PDF equation from five (four 
species and one energy) to three (two species and one 
energy). 

6.8 Revolving Time- Weighted Averaging 

It is noteworthy that although local time-stepping 
seems to overcome some of the problems associated 
with the PDF computations, the application of the 
Monte Carlo method requires the use of a large num- 
ber of particles, because the statistical error associ- 
ated with the Monte Carlo Method is proportional 
to the inverse square root of N m , thereby, making 
the use of the Monte Carlo method computationally 
very time consuming. However, a revolving averaging 
procedure used in Ref. 4 seems to alleviate the need 
for using a large number of stochastic particles, iV m , 
needed in any given single time step. In this averag- 
ing scheme, the solution provided to the CFD solver 
is based on an average of all the particles present over 
the last N av time steps instead of an average solely 
based on the number of particles present in any one 
single time step. This approach seemed to provide 
smooth Monte Carlo solutions to the CFD solver, 
thereby improving the convergence of the coupled 
CFD and Monte Carlo computations. The reason 
for improvement could be attributed to an effective 
increase in the number of stochastic particles from 
N m to N av N m . Here, the solution contained within 
different iterations of the averaging procedure is as- 
sumed to be statistically independent of each other. 

7 DETAILS OF THE COUPLING 
BETWEEN EUPDF-II AND THE OTHER 
SOLVERS 

• The PDF code is designed to be a stand-alone 
computer code which could easily be coupled 
with any of the other unstructured-grid CFD 
solvers (However, some grid-related parameters 
on area vectors, grid connectivity, etc. need to 
be supplied separately). 

• The PDF solver needs information on the mean 
gas velocity, turbulent diffusivity and frequency 
from the CFD solver and the liquid-phase source 
terms from the spray solver. 


• The PDF solver provides the solution for the 
species and energy fields to the CFD and spray 
solvers. 

• It should also be noted that the PDF solver is 
called once at every other specified number of 
CFD iterations. 

• All of the three solvers (LSPRAY, EUPDF-II, 
and CFD) are advanced sequentially in an iter- 
ative manner until a converged solution is ob- 
tained. 

• All three modules (LSPRAY, EUPDF-II, and 
CFD) were coupled and parallelized in such a 
way to achieve maximum efficiency. 

The coupling issues could be better understood 
through the use of a flow chart shown in Fig. 1. 
It shows the overall flow structure of the combined 
CFD, EUPDF-II, and LSPRAY modules. Both the 
PDF and spray codes are loosely coupled with the 
CFD code. The PDF code is designed in such a 
way that only a minimal a minimal amount of ef- 
fort is needed for its coupling with the flow T and spray 
solvers. The present version of the CFD module relies 
entirely on the use of Fortran common blocks for its 
information exchange with other modules. Even this 
reliance should entail only few changes to be made 
within the PDF code for its linkage with different 
solvers. The spray code is also structured along sim- 
ilar coupling principles. 

The flow chart of Fig. 1 contains several blocks 
- some shown in black and/or solid lines and the oth- 
ers in color and/or dashed lines. The ones in solid 
blocks represent the flow chart that is typical of a 
CFD solver. The ones in dashed blocks represent the 
additions arising from the coupling of the PDF and 
spray modules. 

The coupling starts with the calling of the sub- 
routine - spray _pdf_read_parameters, which then 
reads the PDF control parameters from the input 
file, nccjpdf.in of Table 1. The table provides a de- 
tailed description of the following input variables - 
ipdf_mod, ipdLnum, nparti, cmmx, schmitmumber, 
ave.weight, nrollr, pdLcfl, pdf_vardt_fac, pdf_max_dt, 
iran_seed. 

The coupling is then followed by the calling of 
the pdf Jnt _rerun subroutine. It initializes PDF 
computations and, also, it may restart the PDF com- 
putations if needed from the data stored from a pre- 
vious iteration. Similarly, we call spray _int_rerun 
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Fig. 1 The overall flow structure of the combined CFD, LSPRAY, and EUPDF-II solvers. 
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Table 1. ncc.pdf.in file. 

Input file content 

comments 

heading 

title of controlling parameters 

ipdf_mod,ipdf_num 

The variable, ipdf_mod, controls the calls to the Monte Carlo 
PDF solver, eupdf . The PDF solver is called once at every other 
number of CFD iterations as specified by ipdf_mod. 

When eupdf is called, it integrates the PDF computations over 
a fixed number of iterations as given by the number, ipdf_num. 

heading 

title of controlling parameters 

nparti, cmmx, 
schmit .number, 
ave_weight, 
nrollr 

nparti = N m , and it represents is the number of Monte Carlo 
Particles per grid cell. 

cmmx = C 0 , and it represents the constant used in the molecular 
mixing model. 

schmit mumber is used in the calculation of of the 

convection/diffusion step. 

ave.weight has two functions: (1) A positive value for it invokes the 
first averaging scheme and its actual value represents the averaging 
averaging weight (For details, refer to Ref. 5). (2) A negative value 
for it invokes the second averaging scheme described in Section 6 . 8 . 

nrollr = N av , and it represents the number of PDF time steps 
employed in the PDF time-averaging scheme of Section 6 . 8 . 

heading 

title of controlling parameters 

pdf-cfl, 

pdLvardt-fac, 

pdf_max-dt, 

iran_seed 

pdf_cfl represents the constant, Ctf used in Eq. ( 21 ). 

pdf-vardt_fac represents the constant, C% used in Eq. ( 21 ). 

pdLmax_dt is the maximum time step permitted in the PDF 
computations. 

iranjseed is the seed number for use in the random number generation. 
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for the spray computations. It is note- 
worthy that the PDF computations can be 
restarted by reading the data from the restart 
files - nccjpdf -par arris. out, ncc-pdf -result s.db, & 
ncc-pdf -results jave.db. This restart capability is in- 
voked automatically if the needed restart files exist in 
the run-time working directory. Otherwise, the PDF 
computations are initialized to start from the begin- 
ning. Then, the coupling proceeds with the calling 
of the following subroutines: dclr for integrating the 
spray calculations and eupdf for the Monte Carlo 
PDF. The input variable, ipdf_mod of Table 1, con- 
trols the calls to the PDF integration routine, eupdf. 
The PDF solver is called once at every other num- 
ber of CFD iterations as specified by ipdf_mod. And 
every time the subroutine, eupdf, is called it inte- 
grates the pdf computations over a fixed number of 
iterations as given by the input variable, ipdf_num, of 
Table 1. Finally, the coupling ends with the calling of 
a subroutine, spray _pdf .output, which will create 
a set of new restart files. 

8 DETAILS OF THE FORTRAN 
SUBROUTINES 

Table 2 provides a list of all the Fortran sub- 
routines developed as a part of the Monte Carlo PDF 
module. It also describes the purpose of all the indi- 
vidual subroutines. 

9 PARALLELIZATION 

There are several issues associated with the par- 
allelization of both the spray & PDF computations. 
The goal of the parallel implementation is to extract 
maximum parallelism so as to minimize the execu- 
tion time for a given application on a specified num- 
ber of processors [20]. Several types of overhead costs 
are associated with parallel implementation which in- 
clude data dependency, communication, load imbal- 
ance, arithmetic, and memory overheads. The term 
arithmetic overhead is the extra arithmetic opera- 
tions required by the parallel implementation. Mem- 
ory overhead refers to the extra memory needed. Ex- 
cessive memory overhead reduces the size of a prob- 
lem that can be run on a given system and the 
other overheads result in performance degradation 
[20]. Any given application usually consists of sev- 
eral different phases that must be performed in cer- 
tain sequential order. The degree of parallelism and 
data dependencies associated with each of the sub- 
tasks can vary widely [20]. The goal is to achieve 


maximum efficiency with a reasonable programming 
effort [20]. 

In our earlier work, we discussed the parallel 
implementation of a Monte Carlo PDF algorithm de- 
veloped for the structured grid calculations on a Cray 
T3D [6]. The computations were performed in con- 
junction with a CFD solver and LSPRAY. The par- 
allel algorithm made use of the shared memory con- 
structs exclusive to Cray MPP (Massively Parallel 
Processing) Fortran and the computations showed a 
reasonable degree of parallel performance w T hen they 
were performed on a NASA LeRC Cray T3D with 
the number of processors ranging between 8 to 32 [6] . 
Later on, the extension of this method to unstruc- 
tured grids and parallel computing written in For- 
tran 77 with PVM or MPI calls was reported in [7- 
9]. The latest version in Fortran 77/90 offers greater 
computer platform independence. In this section, we 
only highlight some important aspects of paralleliza- 
tion from Refs. [6-9]. 

Both the EUPDF-II and CFD modules are well 
suited for parallel implementation. For the gas-phase 
computations, the domain of computation is simply 
divided into n-Parts of nearly equal size and each part 
is solved by a different processor. Fig. 2 illustrates 
a simple example of the domain decomposition strat- 
egy adopted for the gas-phase computations where 
the total domain is simply divided equally amongst 
the available computer processing elements (PEs). In 
this case, we assumed the number of available PEs to 
be equal to four. However, the spray computations 
are more difficult to parallelize because of its repre- 
sentation based on a Lagrangian formulation. More 
details on its parallelization can be found in Refs. 
[6-9]. 

In the fractional step Monte Carlo method, the 
steps associated with the spatial transport due to con- 
vection and turbulent diffusion as well as the steps as- 
sociated with the transport in the scalar space due to 
molecular mixing, spray contribution, and chemical 
kinetics are advanced in a series of sequential steps 

In the parallel implementation of the Monte 
Carlo PDF code, the inter-processor communications 
are essentially limited to updating the solution at the 
interfaces of the inter-processor domain boundaries 
only once at the beginning of every time step. With 
the domain decomposition adopted, all the stages of 
the fractional step Monte Carlo method, associated 
with such phenomena as the spatial transport asso- 
ciated with convection/diffusion and the transport 
in the scalar space due to molecular mixing, spray, 
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Table 2. Description of EUPDF-II Fortran subroutines. 

Subroutine 

Purpose of the Subroutines 

average (ncyc) 

It calculates the ensemble- averaged temperature and species 
fields based on the averaging scheme selected. 

coeffiv(axyzp,d) 

It calculates the convection/ diffusion coefficients 
as given by c n and c p of Eq. 12. 

convec(axyzp,d) 

It performs the task of moving the particles between the 
neighboring cells during the integration of the convection / 
diffusion step as described in Section 6.2. 

eupdf 

This is the main controlling routine for the Monte Carlo 
PDF solver. When it is called, it updates the PDF solution 
based on Eq. (5) before returning control over to the 
calling routine. 

mimd-pdf 

(nparti,is,phil) 

It facilitates the inter-processor communications by updating 
the solution at the interfaces of the inter-processor domain 
boundaries. 

mimd-pdf_recv 
(i_recvfrom, 
nparti,is, phil) 

It is called by mimd-pdf in order to facilitate the 
inter-processor communications by receiving data 
from the other processors. 

mimd-pdf-send 

(Lsendto, 

nparti,is,phil) 

It is called by mimd.pdf in order to facilitate the 
inter-processor communications by sending data 
to the other processors. 

molmix (freq.dt ) 

It provides the solution for Eq. (22) - the fractional step assoc- 
iated with the molecular mixing as described in Section 6.4. 

outpdfl(ncyc) 

It calculates the residual errors associated with the PDF solution. 

outpdf2(ncyc) 

It creates the restart files for the next continuation. 

pdf_chem2 

It provides the solution based on the integration of Eq. (17) 
- the fractional step associated with the terms arising from 
the chemical reactions. 

pdf_int .rerun 

It initiates the Monte Carlo PDF computations either by 
reading the data from the restart files - nccjpdf -params.out, 
ncc.pdf .results. db, and & ncc.pdf .results .ave.db, or by 
initializing the PDF computations to start from the beginning. 

props-ther 

This routine calculates the following properties for use in the 
CFD solver: (1) temperature, (2) pressure or density depending 
on the flow being compressible or not , (3) specific heat at 
constant pressure, and (4) transport properties. 

ransed(iarg,lu) 

It sets up an appropriate seeding for the random number generation. 

ranu2(xjdim) 

It is a subroutine used in the random number generation. 

spray 

It provides the solution for Eqs. (24)-(25) - the fractional step 
associated with the spray source terms arising from the exchanges 
of mass, momentum, and energy with the liquid-phase. 

vrand(i,x,s,n) 

It generates a vector of pseudo-random numbers uniformly 
distributed on the exclusive interval (0,1) for use 
on 32 bit machines. 
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Fig. 2 An illustration of the parallelization 
strategy employed in the gas flow comp- 
utations . 


and chemical kinetics, lend themselves perfectly to 
parallel computing without any need for interproces- 
sor communications during their integration. There- 
fore, a very high degree of parallelism could easily 
be achieved if enough care is exercised in distribut- 
ing the spatial grid points uniformly amongst all the 
available PEs. Thus, the Monte Carlo simulation is 
ideally suited for parallel computing and the run time 
could be considerably minimized by performing the 
computations on a massively parallel computer. 

9.1 Details of Parallel Performance 

The details of the combined parallel perfor- 
mance of the CFD, EUPDF-II, and LSPRAY codes 
involving several different cases can be found in Refs. 
[6-9]. Here, we only summarize briefly the parallel- 
performance results for two different cases. One is 
a 3D test case and more details on this case can be 
found in the reference [8]. For this case, the calcula- 
tions were performed on a computational grid com- 
prising of 8430 tetrahedral elements and 100 Monte 
Carlo PDF particles per cell. The computations were 
performed on one of the NASA Ames Research Cen- 
ter’s parallel computer platforms called Turing which 


is a SGI Origin work-station with 24 PEs (Processor 
Elements). Table 3 summarizes the CPU times per 
cycle taken by the EUPDF-II, LSPRAY, and CFD 
solvers vs the number of PEs. Both the CFD and 
PDF solvers show good parallel performance with an 
increase in the number of processors but for the spray 
solver it shows reasonable parallel performance. 

Next, we would like to summarize the results 
from [9] showing only the results of the PDF and 
CFD computations. The computations refer to the 
case of a confined swirl-stabilized spray flame. The 
computations were performed on LACE (Lewis Ad- 
vanced Cluster Environment) at NASA LeRC. The 
computations were performed on a grid with a mesh 
size of 3600 quadrilateral elements and a total of 0.36 
million Monte Carlo particles (=100 particles/cell). 
Table 4 summarizes the epu times per cycle taken by 
the PDF and CFD codes versus the number of pro- 
cessors. Both the PDF and CFD solvers showed good 
parallel performance with an increase in the number 
of processors. It takes approximately about 2000 to 
5000 cycles for the computations to reach a converged 
solution. 
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Table 3. CPU time (sec) per cycle versus number of PEs. 



Number of processors 

Solver 

Characteristic 

2 

5 

10 

CFD 

5 steps/cycle 

2.50 

1.25 

0.75 

EUPDF-II 

1 step/ cycle 

6.5 

2.9 

1.9 

LSPRAY 

100 steps/cycle 

1.70 

0.64 

0.53 


Table 4. Cpu time (sec) per cycle versus number of PEs on LACE Cluster. 



Number of processors 

Solver 

Characteristic 

2 

4 

8 

16 

EUPDF-II 

1 step/cycle 

2.30 

1.35 

0.75 

0.44 

CFD 

5 steps/ cycle 

3.55 

1.90 

1.10 

0.60 


10 A SUMMARY OF SOME RECENT 
VALIDATION CASES 

The following eight cases were validated: 

1. A supersonic axisymmetric jet involving hydro- 
gen/air combustion. 

2. A coaxial supersonic burner involving hydro- 
gen/air combustion. 

3. A planar high-subsonic reacting shear layer in- 
volving hydrogen/air combustion. 

4. A reacting methanol spray with no-swirl. 

5. A non-reacting methanol spray with no-swirl. 

6. A confined swirl-stabilized n-heptane reacting 
spray. 

7. An unconfined swirl-stabilized n-heptane react- 
ing spray. 

8. A confined swirl-stabilized kerosene reacting 
spray. 


The experimental data for the first case was pro- 
vided by Evans et al [21] and for the second case it 
was provided by Cheng et al [22] . The validation for 
Cases 1 and 2 was described by Hsu et al in [4] and [5] 
and the comparisons showed reasonable agreement. 
And the experiment ah data for Case 3 was provided 
by Chang et al [23] and its validation involving a pla- 
nar high-subsonic reacting shear layer was reported 
by Liu & Raju [24]. 

The experimental data for Cases 4 & 5 was pro- 
vided by McDonell & Samuelsen from the University 
of California at Irvine [25]. Both the cases are with- 
out swirl; one is a reacting case and the other is non- 
reacting. The data for Cases 6 & 7 was provided 
by Bulzan from the NASA Glenn Research Center 
[26-27]. Both the cases are swirl-stabilized reacting 
cases, one is an unconfined flame and the other is 
confined. The data for the last case was provided by 
El Banhawy & Whitelaw from Imperial College [28]. 
It is a confined swirl-stabilized kerosene spray flame. 
Here, we would like to provide a brief summary of the 
validation cases but a detailed presentation of the re- 
sults and discussion can be found elsewhere in the 
papers [6-9]. The comparisons involved both gas and 
drop velocities, drop size distributions, drop spread- 
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ing rates, and gas temperatures. The results were in 
reasonable agreement with the available experimental 
data. The comparisons also involved the results ob- 
tained from the use of the Monte Carlo PDF method 
as well as those obtained from a conventional CFD 
solution without the Monte Carlo PDF method. For 
the case of McDonell & Samuelsen’s reacting spray 
flame, the detailed comparisons clearly highlighted 
the importance of chemistry /turbulence interactions 
in the modeling of reacting sprays [8]. The results 
from the PDF and non-PDF methods were found to 
be markedly different with the PDF solution provid- 
ing a better approximation to the reported experi- 
mental data. The PDF solution showed that most 
of the combustion occurred in two distinct flame re- 
gions with most of the combustion occurring in a 
regime that is mostly diffusion controlled and the sec- 
ond regime showing the characteristics of a premixed 
flame. However, the non-PDF predictions showed in- 
correctly that most of the combustion occurred in a 
predominantly vaporization-controlled regime. The 
Monte Carlo temperature distribution showed that 
the functional form of the PDF for the temperature 
fluctuations varied substantially from point to point. 
The results brought to the fore some of the deficien- 
cies associated with the use of assumed-shape PDF 
methods in spray computations. 

11 CONCLUDING REMARKS 

• This manual provides a complete description 
of EUPDF-II - An Eulerian Monte Carlo PDF 
solver developed for application with parallel 
computing and unstructured grids. 

• We have extended the joint scalar Monte 
Carlo PDF method to two-phase flows and, 
thereby, demonstrating the importance of chem- 
istry/turbulence interactions in the modeling of 
reacting sprays. 

• The method outlines several techniques designed 
to overcome some of the high computer time and 
storage limitations associated with the Monte 
Carlo simulation of practical combustor flows. 

• It provides the user with a basic understanding 
of the PDF formulation and the EUPDF-II code 
structure, and complete details on how to couple 
the PDF code to other CFD codes. 

• The basic structure adopted for the grid repre- 
sentation and parallelization follows the guide- 
lines established for NCC. 


• Based on the validation studies involving several 
spray flames, the results were found to be en- 
couraging in terms of their ability to capture the 
overall structure of a spray flame. 

• The source code of EUPDF-II will be available 
w r ith NCC as a complete package. 
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